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FOREWORD 


The present report describes a rocket nozzle turbulent, boundary layer 
heat transfer rate parametric study using three different computer programs. 

This investigation, entitled ROCKET ENGINE THRUST CHAMBER HEAT TRANSFER 
CALCULATION AND ANALYSIS, was conducted for the Propulsion and Thermodynamics 
Division of MSFC, National Aeronautical and Space Administration under the 
grant No. NGR01-0Q1-024 with Klaus W. Gross as technical monitor. 

Hrishikesh Saha of Alabama Agricultural and Mechanical University was the 
principal investigator. 

This research program contributed extensibely to improve faculty and 
student-research capability at the Alabama A§M University and this support 
by NASA-MSFC is greatly appreciated. 

The author wishes to acknowledge the helpful assistance and advice 
received from Mr. Klaus W. Gross' and Alfred Krebsbach, NASA-MSFC. The 
author would also like to express his appreciation to Mr. Marion I. Kent, 
Assistant Director, University Affairs, NASA, MSFC, for 'his effort in 
making this grant possible. 
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ABSTRACT 


A parametric study of the heat transfer rate along the wall of a 
rocket nozzle is presented. The influences of different parameters; ' 
Laminar and turbulent Lewis number, mixture ratio, initial wall temperature 
distribution, and eddy viscosity, are considered in this study. 

The numerical evaluation of these influences on heat transfer rate 
is done by using three different compressible, reacting Laminar and 
turbulent boundary layer computer programs; MABL (Mass Addition Boundary 
Layer Program) , MABL-KI (MABL program is modified to include turbulent 
kinetic energy equation) , and BLIMP (Boundary Layer Integral Matrix 
Procedure) . 

This study also provides ah' excellent opportunity to evaluate the 
efficiencies of these three computer programs and to suggest one of them 
for future computational purposes. As described in the recommendations, 
the suggested program needs further improvements, in its mathematical 
model and numerical solution method, to be used as a standard computing 
procedure for comparison with future experimental results. 
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SECTION I > - 

INTRODUCTION 

Boundary layer behavior along the walls of a rocket 
nozzle plays an important role in the performance of a 
rocket engine. The shear layer determines part of the 
thrust loss of the nozzle and the energy layer controls 
the heat transfer to the wall and the wall temperature 0 
There is a need for a computer program which can calculate 
boundary layer properties and effects for flows with large 
pressure gradient, chemical reactions, and a wide variety 
of wall conditions o 

Accurate heat transfer prediction for regeneratively 
cooled thrust chamber is of vital importance, since the 
increased energy level of the coolant has a direct impact 
on the engine performance* Heat transfer rate has also an 
effect on the contour design and material selection of the 
thrust chamber wall. 

i 

The primary purpose of this report is to present the 
results of a comparative study of boundary layer heat trans- 
fer rate at a nozzle wall, computed by the computer programs 
presently available at the MSEC Propulsion and Thermodynamics 
Division and also to select one of the programs suitable for 
further studies of the SSME characteristic s 0 

The computer programs presently available at NASA-HSFC 
for prediction of heat transfer rate in a rocket nozzle bound- 

■i 

ary layer are 
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a) MABL ( Mass Addition Boundary Layer, Ref. 11&X2). 

b) MABL -KjE ( MABL modification to include turbulent kxnetic 

energy equation., Ref* c 19 ) * 

c) BLIMP ( Boundary Layer Integral Matrix Procedure, Ref .23-25 ) . 

Short descriptions of these programs with respect to their, 

- basic assumptions made to set up the governing equations 
of the boundary layer flow including discussions of the con- 
servation equations, turbulence model, coordinate transforma- 
tion, and boundary conditions, and 

- numerical techniques and analysis which led to the computer 
solutions, 

are given in section 2o 

For the sake of completeness, a short derivation of the 
steady state compressible turbulent boundary layer aquations 
for two dimensional and axisymmetric flows of chemically re- 
acting mixtures is given in the appendix. The turbulent kin- 
etic energy equation to compute eddy viscosity is also includ- 

ed in the appendix « 

Section 3 contains the various turbulence model, definitions, 
and presently available computational formulas of the laminar 
and turbulent transport coefficients and parameters; coeffic- 
ients of molecular, -viscosity , jl , thermal conductivity, K, 
and diffusion, 3>t j , and eddy viscosity £, eddy conductivity 
Kp, eddy diffusivity Df, and Prandtl number, P„ Lewis 'number, 

L, and Schmidt number, S. This section also contains a deri- 
vation of the heat flux equation used in different computer 


programs * 
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A parametric study of the heat transfer rata along the 
wall of a rocket angina thrust chamber using different com- 
puter programs is given in section The effects on the 
heat flux due to the variation of the parameters; Lewis 
number. Mixture ratio, initial wall temperature distribution, 
and eddy viscosity are shown in the figures 0 

Section 5 presents a comparative study of the available 
computer programs, MABL, MABL-XS, and BLIMP, and suggests one 
program for future computation with several modifications and 
improvements 0 
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SECTION II ' 

BOUNDARY LAYER COMPUTER PROGRAMS 

The computer programs, presently used at MSFC to compute turbulent boundary 
layer properties in a rocket thrust chamber, are described below: 

MABL (Mass Addition Boundary Layer) Computer Program . 

This effort was carried out (Ref 11 § 12 } for the purpose of developing 
a practical and immediately useful state-of-the-art engineering method for 
analitically predicting the boundary layer performance losses in liquid 
rocket engines; in particular for the space shuttle engine configurations 
currently under consideration by NASA. Therefore, the analysis and 
resulting computer program have been oriented towards the analysis of 
high pressure hydrogen- oxygen engines. 

The boundary layer equations for compressible turbulent flow, 
derived from the time dependent Navier-Stokes equations using the Reynolds 
time-averaging procedure and the usual boundary layer order of magnitude 
assumptions in a curvilinear coordinate system, were considered for solution. 

The analysis and program use equilibrium chemistry since the flow in 
the high pressure H o -0^ engines will most probably be very close to the 
equilibrium. The computer program is set up and dimensioned, at present, 
to handle the hydrogen- oxygen system and considers two elements (H,0) and six 
species (H^, H, OH, 0 ? , 0) . The equilibrium state calculations in the 

present analysis are performed by selected subroutines from the One-Dimensional 
Equilibrium (ODE) JANNAF Performance Program; given the element mass fractions 
and two thermodynamic state variables this program solves for the species 
mass fractions and all the other thermodynamics state variables. 
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As a result of the parabolic nature of the boundary layer equations, 
initial profiles are required to begin the solution. 

The numerical method chosen for the solution of the boundary layer 
equations is a Crank-Nicolson implicit finite difference technique. In 
applying this finite difference technique the difference analogs of the 
equations of motion are linearized and uncoupled and solved using a 
tridiagonal matrix inversion algorithom. 

In formulating the turbulent equations of motion it was assumed that 
the turbulent flux terms could be related to mean, time averaged 
quantities, through the use of the phenomenological mixing length-eddy 
viscosity concept. The turbulence eddy viscosity model, used here, was 
developed and extended by Cebeci (Refs. 17, 18). This model uses a two- 
layer representation of the eddy viscosity. In the inner region, closer 
to the wall, the eddy viscosity is based on Prandtl’s mixing length theory, 
as modified by Van Driest to account for the damping effect of the wall, and 
as extended by Cebeci to include wall mass transfer, compresibility and 
pressure gradient effects. In the outer, wake-like, portion of the boundary 
layer, Clauser's form of the eddy viscosity, modified to include an 
intermittency factor, is used. Turbulent Prandtl number formulation 
developed by Cebeci (Ref. 17) is used here. The equations and computer 
program have been kept general to allow for an arbitrary variation of 
turbulent Lewis number. No attempt has been made to calculate or model 
Lewis number. 

MABL-KE (MABL with Turbulent Kinetic Energy Equation) Computer Program . 

The MABL program described earlier was modified to include turbulent 
kinetic energy equation for solving the compressible turbulent boundary 
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layers along the wall of a rocket thurst chamber (Ref. 19). This inclusion 
of the turbulent K. E. equation provides better results than the method 
utilizing the Van Driest -Clauser hypothesis; this also the character of 
chemically reactive flow with mass injection and can exhibit the past his- 
tory of boundary layers. Results show the variation of eddy viscosity in 
flow direction including the familiarization tendency in the nozzle conver- 
gent section. The results appear to be in good agreement with available 
experimental data. The eddy viscosity decreases significantly in the re- 
active hydrogen-oxygen boundary layer when a small amount of hydrogen is 
injected from the wall; 

The equations and boundary conditions for the compressible turbu- 
lent boundary layer are included in the appendix for completeness sake. 

The numerical solution technique is similar to that of MAHL as 
described earlier. 

/ 

BL I^_ (jounda,ry Layer Integ ral M atri x Procedure) Computer Program . 

ELIKP version J serves a.s a boundary layer prediction computer program 
for rocket nozzle flow. The references 13, 25, 29, 30, 31, 32 describe a 
mathematical model and numerical solution of nonsimilar laminar or turbu- 
lent multicomponent chemically reacting boundary layer flows with arbitary 
equilbrium chemical systems, unequal diffusion and thermal diffusion co- 
efficient for all species, and a variety of surface boundary conditions. 

±ne equations of motion are developed for axisynetric or two-dimen- 
sional turbulent flow including transverse curvature effects. The usual 
turbulent flow technique of breaking the species, velocity. 
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and enthalpy fields into mean and fluctuating components,, time averaging, 
and making appropriate order of magnitude approximations results in the 
governing equations . The turbulent transport terms are expressed in the 
Boussines^ form, i.e. eddy viscosity, eddy diffusivity, and eddy 
conductivity. Hence, all the terms in the equations are time- averaged 
quantities. In the order-of-magnitude arguments, terms of the following 
types have been eliminated; 1) triple correlations, 2) derivatives of 
turbulent correlations parallel to the wall, and 3) correlations involving 
turbulent components of molecular transport mechanisms. Also necessary in 
the mathematical formulation of the problem is the specification of the 
molecular transport properties, equation of state and equilibrium relations 
for the multicomponent gas, and a descritpion of the eddy viscosity, 
conductivity, and diffusivity. 

The boundary layer fl : ow is divided into a wall region and a wake 
region. A mixing length formulation for turbulent shear stress was chosen 
for the wall region. In the wake portion of the boundary layer, a constant 
eddy viscosity model (Clauser's equilibrium boundary layer expression 
(Ref. 33)) is used. 

A bifurcation approximation to binary diffusion coefficients 
introduced in (Ref. 20) utilized herein permits expilicit solution of the 
Stefan-Maxwell relations for the diffusion flux in terms of gradients 
and properties of species and of the system as a whole. 

A modified Levy-Lees coordinate transformation is used to reduce the 

partial -differential equations of the boundary layer to total-differential 
equations . 



8 


The solution procedure uses splined gradratics or cubics to represent 
velocity, enthalpy, and species concentrations between nodes which give 
smooth variations of the dependent variables and their derivatives through 
the boundary layer, therefore allowing significantly fewer -nodes for a 
given accuracy. This advantage is particularly valuable for turbulent 
boundary layers, where the very large gradients in' the surface normal 
direction would require a great number of nodal points for an accurate 
solution with typical finite difference representations. In the streamwise 
direction, derivatives are expressed by ordinary implicit difference 
relations using linear or quadratic curvefits of the dependent variables. 
The resulting set of algebraic linear and nonlinear simultaneous equations 
is solved iteratively by the general Newton- Raphson technique. 
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T KA.H 3? 0 NT ? NO PE NT IE 3 


histories lly there have been two different approaches to 
developing the equations of gas dynamics, i.e. equations for 
the conservation of species, conservation of mass, conserva- 
tion of momentum, and conservation of energy along with the 
equation of state for the gas mixture. 

One method is the phenomenological approach wherein 
certain relations between shear and strain, heat flux and 
temperature gradient, and diffusion flux and concentration 
gradient are postulated and the equations then developed using 
the laws of classical mechanics and heat flow* The method 
leaves the transport coefficient s, that is, the constant of 
proportionality between shear and strain, heat flux and temp- 
erature gradient, and diffusion flux and concentration grad- 
ient, undefined 'with no method other than direct measurements 
available to determine their values. These transport coeffi- 
cients are identified by the symbols^, K, and Dfcj for laminar 
flows, where, for example (Ref. 3 ), 

viscosity 

~ Coefficient of thermal conductivity 


and 


. Vue-, _ 

‘J '■HC-./Zi ~ “• 


nary diffusion coefficient 
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where - shear stress on surface perpendicular to y 

azis acting in x direction in cartesian co- 
ordinate system, 

# 

~ heat flux in y direction of cartesian coor- 
dinate system, 

' “ component _ of diffusion velocity due to con- 

centration gradients in y direction of carte- 
sian coordinate system for a species i diffus- 
ing into a mixture of species i and j, 

c 

I ~ ?i/? 9 mass fraction of species i, 

T and S are the temperature and density of the gas mix- 
ture, respectively. 

The alternative approach is that of the mathematical theory 
for nonuniform gases, essentially a kinetic theory approach (Ref, 
21), This method yields the fluid-dynamic equations viith the 


transport coefficients defined in terms of certain integral 
relations which involve the dynamics of colliding particles., 
Some empiricism is involved in specifying the interparticle 
forces needed to evaluate the collision integrals but even 
hero, in principle, recourse to theory will narrow the margin 
of uncertainty involved. The transport coefficients jri, K, and 


D{j arc found to be funtions of the 
gas -species molecular weights, and c 


gas -mixture t ©mperature , 
ertain parameters of the 


1 nt e rp art ids f o r c e f i eld. 

The conservation equations can often be written in forms 
involving d i men s i unless ratios or various t ransportcoeff iciert s. » 
These dimensionless transport parameters arc defined as follows: 
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The PR AN .DTE 




union i 


r-ouxh measure of the relative 


importance of momentum transfer and 
P r -juC F j/}< pfnere C pj; = Ci Cp L 
h Specific heat 

r t j - 


heat transfer, is defined as 
= frozen Specific heat 


It can be shown from the- energy equation (Hef. 34) that 
K’ (^T/fh) = Conduction _ _JL -where 
J*i It {(cjllfdy} Shear ‘work P r B 

E -Li e/he is called Eckert’s numbers 
For a pure polyatomic gas Prandtl number is given by the 
Euekon form, P== 4 T / (9 T -5)« From this equation it is seen 
that ? r , which is approximately 3/4 for air ( V =1 ,i|- ) , varies 
from 2/3 for monatomic gases { t =5/3 ) to unity as In 

some systems ( e ♦ g * , liquids) the Eucken , formula is not valid 


and P r c.an differ considerably from unity. 


The 8 C Ill'll FT NUMBER, a rough measure of the relative importance 
of momentum transfer and mass transfer, is defined for a b in ary 


mixture by, 

8 =y } i /SD {2? from Energy equations ( Ref, 340 it can be 
shown that § hi(<dCi fy y ) = diffusion -^ 1 

J*L LL ( *cs ui®y ) , shear work 3E 

If S Is large, shear work predominates over diffusion. In 


multicomponent systems, Schmidt numbers may be defined for 
each pair of species * As in the case of Prandtl numbers, 
somewhat less than unity. 


is often 
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Ihe a-Frid is trio ratio of the energy transported toy 

conduction to that transported by diffusion, and is defined 
for a binary mixture by 

1 ” ^ jjf £ Q&£ ~ Br ^ di ffusion - ^ k- hz(fo C:/0N ) 

K S ^ conduction K ( qyy /'&y } 


As with the Schmidt number, Lewis numbers may be defined for 
each pair of species in multicomponent mixtures. 

In many systems, L is very nearly unity; it is often 
slightly greater than unity in combustible gas mixtures. 

The approximation L= 1 is frequently very helpful in theoret- 
ical combustion analysis. 


LAMINAR TRANSPORT PROPERTIES 

The ideal gas calculations performed to date wore basically 
for program check-out purposes and restricted to a one com- 
ponent gas so the Lewis number is not required. The Prandtl 
number is assumed to be a constant. Its value can be selected 
to approximate the gas boing considered. Most of the ideal gas 
calculations have been for air and currently the viscosity 
calculation is oasecl on Sutherland 1 s law expressed as 

jvl = 2.27 ‘ lo' 8 . T 1/2 / ( 1 * 198,6/T) 

For gases other than air a different viscosity formulation 
s ho u 1 d b e u a & d o 

'.’..0.6 calculation of the laminar transport properties for 
a mixture or gases (ror example : H? -Op system } is based on 
the following: 

The viscosity of the mixture is calculated from Tilke T s 
s o mi - eiap eric, a 1 f o rmu 1 a ( Re f . 20). 

ORIGINAL PAGE IS 
#00B QUAUTYi 
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Vs 


n< 


j L = 


L~ I 


«r-' 


‘J * 


- 1 

J : ) 1 


where h g is the number of species, Xj_ is the raol fraction of 

species i; that -is, X^ = C 1 / ( * 21 Cl / Hi ), M i is the 

i J 

viscosity of each species 1; 

J*i = 266.93 X 10“7 (Hi T)' 5 / ( 0 2 .Q_i2,2 ) 

here = molecular weight of species i v 
= collision diameter, 

T = temperature 

-Oi 2 # ^- ) * v - > 2 ^ / { _qJ 2 , 2 ) ) rigid sphere 

_qJ 2 ,2) - collision integral ( Ref. 21) 

snd Xlj =1} + (J H 2 /{(|i 3/2 ( 1 + Ml / Mj 
Por simpler and faster computational purposes the values of 
J * i are taken from reference U 4 . 0 ] and stored in tabular form. 

The more accurate predictions of Surehla foreland P r are given 
in reference (1+2). 

THERMAL CONDUCTIVITY. Por a monatomic gas of species i the 
Chapman -Enskog theory yields ( Ref. 21) 

% - ii hi ■ 

4 Ki 

•where K^~ thermal conductivity 

R ~ universal gas constant 

Mason and Saxena". ( Ref. I4.I) give the thermal conductivity of 


a mixture of polyatomic gas by 
v$ 

ir _ r If ■ . / . . -,065 

J = < 
0 >» 


r K; L> + ' 

t-; 


u ,y 

tw^xj/xJ 


-1 


where ZCi - l£ Jt ± ( „35J>. C pi + .115 • R / Kj. ) 

4- 
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DIFFUSION COEFFICIENT, 
coefficient is 

D 12 = D 21 = x 

where S' = + S^) 


The expression for the binary diffusion 



2 & - Transport property collision dia. 


p = pressure 
T = temperature 


The effects of multicomponent diffusion can be considered by 
using a bifurcation approximation as given in reference [26]. 

The frozen specific heat of the mixture is 

c pi - £ 

and the laminar Prandtl number can then be calculated as 

Pr = f- C M / K 

and similarly the laminar Lewis number for a binary mixture may 
be computed as 

L = SVuCtf/K 


TURBULENT TRANSPORT PROPERTIES 

The turbulent equations of a boundary layer contain terms 

which represent the transport of mass, momentum, and heat due 

to the turbulent fluctuations and are identified as follows: 

" C § V) Q " — ’ turbulent mass transfer 
~ ~7 f 

“ C §v U ~ turbulent momentum transfer 

7— T 

- C S v/J n l — ' tu rbu 1 q nt heat t ran s f e r 
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The eddy viscosity, 8 , eddy conductivity, ICp, and eddy 
diffusivifcy, Dv, are defined by assuming that the turbulent 
flux terms can be related to mean, time averaged quantities, 
through the use of 'the phenomenological mixing length-eddy 
v ;I s c o s i t y c o nc e pt a. s fo llo v; s : 

ecs.P = - 


c s v)' u' 


0u/^' 


i< (s ; p - _ infh 

T ® T 

g pT (s,q - - Y gjfy' 

<0C- fo* 

L 

Furthermore , turbulent Lewis number and Prandtl number 


are defined, using the above, as 


ijiji k 


, c Pt g d t 


K-, 


and Prp = 


Cp; £ 

Kt 


SDDY VISCOSITY. 

The details of CobeciJs extended eddy viscosity model is 
given in reference (l8)» The model uses a two -layer represent- 
ation of the eddy viscosity. In the inner region, closer to 
the -wall, the eddy viscosity is based on Prandtl 1 s mixing length 
theory, as modified by Van Driest to account for the damping 
effect of the wall, and as extended by Cebeci to include wall 
mass transfer, compressibility and pressure gradient effects. 
Thus, the eddy viscosity in the inner region is given by 

rr g L~ | | . where the mixing length, 1, is 



exp (~y/A) ] 
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the Van Driest 1 s "damping factor". A, is defined as. 

A= [26 ft]/ [ (T w § ft n] 

where = shear stress at the wall 

and the factor, N, which accounts for pressure gradient and 
mass transfer effects is given by 


N 2 


dP Jt i 
d X 




exp 


H.8C*Y)wJ ll w 


-f- exp 


C§v; w c M l . ' CZwVw )' 12 ./ 1 
11.8 C^J w J^w 


-J 


cr„<rj "’jn 


If (• § V ) equals 
<v 


zero the above equation becomes a degen- 
erate form and N must be calculated as 


h 2 = lni.BdP I's to -m. ( 3 >is 


In the outer, wake -like, portion of the boundary . layer , 
Clauser * s form of the eddy viscosity, modified to include 
an intermittency factor, is used. The outer eddy viscosity 
is given by 


£ c - 0*0168 


u e * 5.5 ( y /§ ) 6 ] 


-1 


where the term in brackets is an approximation to Xlebanoff » s 
error function intermittency relationship and _ is the in- 

L T3 C 

c o rap res sib 1 e d i s p 1 a c e me n t t hie kn ess d e fin ad as 


&I* e =J 0-^Me)dy 

O 


The eddy viscosity for the inner region is used from the 


rail ou 


t ward until the height at which £ ~ f ■ 1 

- o 


is reacned. 


roi 


that point, -to the boundary layer edge, the outer expression 


n 
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for eddy viscosity is utilised. 

For an external air flow the results using Van. Driest - 
Clausen eddy viscosity model are in agreement with experi- 
mental measurements ( Ref. 18 ) . However, for an internal 
flow as in a rocket thrust chamber with a highly cooled wall, 
C-ebeci's modified model does not show the laminari nation 
tendency in the nozzle convergent section and the increase 
of eddy viscosity in the divergent section of the nozzle 
even for air flow (Ref. 19). In order to include the effect 
of the past history of flow and the strong temperature var- 
iation across the boundary layer due to chemical reaction 
' . • 
or wall cooling, the turbulent . kinetic energy equation has 

been introduced for the solution of the Reynolds stress (Ref. 
19, also see Appendix). This method obtains eddy viscosity, 
£ 0 , in the wake region, and uses the formula introduced 
earlier for the inner eddy viscosity, <fh . 


TURBULENT PRANDTL NUMBER . 

The experimental data (Ref, 16) have shown that the turbulent 
Prandtl number varies considerably across the boundary layer. 
To account for this type of variation the 'turbulent Prandtl number ■ 
is computed using Gebeci’s formulation (Ref. 17) as 

P T = K-n. [l ~ exp { / £% p-e.tp(-y+|y2 /3 + )]j 
•where lh, - .1, ,44 are the momentum and enthalpy Prandtl 




length constants, respectively, and 


j * ~ 7 Li# / 

where V - Hi nematic vis go sit v 


6 ( § u /§) 2 ( 1 


^JSSg , 



IS 


B + 


3k ( %/ <S P (ft/ jy l w ) i 


u/\. _ 


( Tr w /f„) * 


at the wall, y=0, the above equation reduces to 
P T = fi&n B + j / [K n A+ p y j 


Hot too much is known about the turbulent Lewis number* and 
eddy conductivity, K<p, and eddy diffusivity. Dp, at this time. 

In the calculations L<p has been assumed constant close to unity. 
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HEAT FLUX tj 

The local heat flow per unit time and area (heat flux) in a pure 
substance is given by the Fourier* s law of heat conduction: 

■q « -k VT (I) 

It states that the heat flux vector q is proportional to the tem- 


perature gradient VT and is oppositely directed; k is the proportionality 
constant or thermal conductivity. Thus in tan isotropic medium heat 
flows by conduction in the direction of steepest temperature descent. 

In a moving fluid q represents the flux of thermal energy relative to 
the local fluid veclocity* 

For mixtures there are, in addition to the conductive flux, con- 
tributions resulting from the interdiffusion of the various species 
present and the diffusion-thermo (Bufour) effect. The total energy 
flux relative to the mass average velocity is then given by 

q -q Cc)^(d) + -£ (x) (2) 

(The radiant energy flux may be handled separately and will be 

described later). Here in equation (2) q£ c )=-k VT is the conductive 
energy flux, as defined in equation (1), and k is the instantaneous 
local thermal conductivity of the mixture. The energy flux caused 

by interdiffusion is defined for a fluid containing n species by the 








J i 


expression: 


C3) 
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Here H i is the partial molal enthalpy of the it hr species, and 
is the diffusion ijiass flux in a multicomponent system# The Bufour 
energy flux q is quite complex in nature and is usually of minor 
importance# The explicit form of the Dufour-ef feet term in multicomponent 
gas mixtures has been discussed in [Ref. 2l], 


In order to compute heat transfer rate in a reacting turbulent 
boundary layer using the computer program MABL (Mass Addition Boundary 
Layer) [Ref. 11&12], the approximate heat flux relation given in this 
program is derived below; "The peculiar veclocity of species 1 is defined 
as 


_ _ 

v* = i>. - It 

* l 0 

and is the veclocity of a particle of kind i with respect to a coordinate 

system moving with the mass average veclocity of the fluid, where u . 

i 

is the random veclocity of the i th. species*. 

The diffusion veclocity of species i with an average veclocity "£?• 
in a mixture flowing with a mass average velocity u* 0 is defined as 

\ = ^r ^ 0 m 

This veclocity is zero when only one kind of particle is present# The rm*9$ 
flux (per unit time and area) of species i across a surface moving with 
the mass average veclocity of the mixture is given by 

Mass Flux^ = v* , where density of species i 

Then the enthalpy flux across this surface due to the mass flux of 

* 

specieslis 

Enthalpy f lux . - §^V. K< ($ ) 

Where h* * enthalpy per i species# 
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The total enthalpy transport '( of all kinds of species) per unit 
time and area is 

Z h. (63 

since mass fraction C.= j f 

£ *1% c i b i (?) 

2*1 t=l 

Fick* s law gives the relation (approximate formula) 

\> = - D VClog C* ) 

Where D = diffusion coefficient or, 

Vf “ DVC i (S) 

Substituting relation (8) in relation (7) the energy flux caused by 


interdi f fusion 


of species in a gas mixture is derived, i*e., 


q Cd) = - §D?h;VC; 

The total heat flux^tothe mass average veclocity as in equation 
(2) neglecting tha.Bufour effect is given by 

t = q (c) + q (d) 

q=-K7T - ?D Eh: 7C. (.10) 

U/ *• 1 

To express equation (10) in terms of enthalpy^T may be transformed as 

n 

Since, total enthalpy h — Z C. h- (ll) 

i* i 1 1 

grad h = Vh = (C-Vh ; * h.FC ; 

oL ' *t o 

Now, ’7h i = = 3(/CpjdT+ h; ) ?T 


v VT 


where the relation h* — fc„_>dTd- h, was used 

rl P 1 t 


= Specific heat & h, = heat of formation of species i 



Substituting Egn (13) into Egn. (12) we get 

C + h. VC* ) 
t~! ^ n 

= C V T + £ h; VC* j where C = 

Pi i,| 1 1 P$ 

or, VT = i (7h - £ h^VCi ) l )U ) 

Substituting Egn (14) in Egn, (10) and rearranging we g.et 

"5 = -^ (7h-|h;VC; ) -fD^h/VC; 

Pi t*l ‘ i*/ 

== - y ~ [Vh +?h.VCj (?JL£ej- l) ] 

/* c pf t=j 1 K 

t= -iL [Vh+ (l, - l) j?h; VC; ] (<5> 

p r i*\ c L 

where Prandfcl number P r = Cp y-l and Lewis number L = g D C p 

T R. 

This is the heat flux equation used in the. computer program KABL. 


M3 
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SECTION IV 

RESULTS AND DISCUSSIONS OF THE PARAMETRIC STUDY 

The computer programs, MABL, MABL-KE, and BLIMP are used to predict 
the heat transfer rate at the walls of NASA's Space Shuttle Main Engine 
(SSME) and RL-10 Engine thrust chambers. The results are shown in figures 
1 through 7 for various parameter values. 

The SSME characteristic parameters are: 

Nozzle geometry - area ratio, E=77.5, throat radius, R =0.4294 (ft). 
Stagnation conditions for normal power load- 
temperature, T o =6000° (R) ’’ 
pressure, P q =3000 ( psia) 
mixture ratio, MR=6.0 for O 2 /H 2 fuel system. 
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The PL-10 Engine characteristic, parameters are: 

Nozzle geometry - area ratio, 3-o7 

throat radius, Rt-0.2lUl6(ft) 

Stagnation conditions 

temperature, T =£000° (R) 
o 

Pressure, P =1*00 (psia) 
o 

Misture ratio, KR*J. 0 for °2/ h 2 fuel system. 

The heat transfer rate predictions at the wall of the SSME thrust cham- 
bers using KABL, MAEL-KE and BLIMP are shown in figurel, where the abscissa 
represents the wall contour length. The thrust chamber radius variation is 
shown against the contour length to aid the visual apprehension of the heat 
transier rate predictions along the wall contour: The predictions have the 
maximum heat transfer rate slightly upstream of the throat as was found by 
a few experimental results. . Since ''the heat flux computed by the HABL pro- 
gram differs slightly from the results of MABL-KS, only the computed values 
near the throat region are shown in the figure 1. 

EFFECTS OF LEWIS NUMBER ON HEAT FLUX 

To determine the influence of Lewis number on heat flux the computer 
program MABL-KE was used for SSME & KL-10 engine configurations, since the 
BLUiP program does not explicitly express its governing equations in terms 
of lewis number. Since the MARL-ICE does not compute Lewis number, fixed 
values about 1.0 as input to the program were chosen. Figure 2 represents 
the heat flux variations for L-l . 0=L^ for a fixed mixture ratio of, . 

Even thus small change in Lewis number causes visible difference in the heat 
flux 
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distributions. This influence is more prominent at the throat region, which 
was expected because of the sharp variations of the flow parameters there. 

Heat flux increases with Lewis number for a fixed mixture ratio. KABL-KE 
faces computational difficulties for higher Lewis numbers . 

The BLIMP program considers multi-component diffusion effects by using 
a bifurcation approximation, which means that the Lewis number, considered but 
not explicitly computed by BLIMP, would be more accurate for a boundary layer ( 
flow analysis than the fixed Lewis number as assumed in KABL & KABL-KE. 

EFFECTS OF MIXTURE RATIO OH HEAT FLUX 
'The influence of mixture ratio on heat flux was determined by using 
BLIMP for the SSMB configuration. Figures 3,k, and 5 show the influences on 
heat flux due to the mixture ratios, MR^U, 5, & 6,U£, each for a fixed wall 
temperature distribution respectively. In these figures heat flux is plotted 
against the axial thrust chamber distance for the SSHE configuration only. 

The heat flux at the wall increases with the increasing mixture ratio and is 
more prominent neon the throat region, as expected. Figure h . shows the 
nominal wall temperature distribution which is a required input for computation. 
This temperature profile has been predicted analytically by a different computer 
program utilizing test results from other similar thrust chambers. 

EFFECTS OF INITIAX WALL TEMPERATURE DISTRIBUTION ON HEAT FLU X ■ 

Since the initial wall temperature distribution must be given to initiate 
computation and since this temperature profile is a predicted quantity con- 
taining same uncentainty, it is of great importance to observe its influence 
on the heat flux. Figures 6 and 7 reproduces . • 
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the heat flux results computed by using the BLIMP program for the S3MS thrust 
chamoer. Three different wall temperature profiles, T wall, (T wall + 200 °R), 
ana (T vs 11 - 20O°E) were used to compute the heat flux for several fixed mix- 
ture ratios respectively. Heat flux increases with decreasing initial wall tern- 

psrature and this influence reaches its mszujrmm near the throat region as ex- 
pected . 

AFFECTS OF SDDr VISCOSITY KODSL AND PRANDTL NUMBER ON HEAT FLUX 

- Z visgo sity ! Hash of the computer programs, MABL, KABL-KE, and BLIMP, 
uses different eddy viscosity models. KABL uses the Van Briest-Clauser 
formula .modified by Cebeci (Ref. 1?). KAHL-KS modifies the MABL formulation by 
introducing the turbulent kinetic energy equation to find eddy viscosity in 
the -wake region (Ref. -19). 

BLDIP uses the Qlaussr formula without any modification to compute eddy 
viscosity, Omori et al., in reference (1?) compares eddy viscosity profiles 
across boundary layers, computed by using MAHL-XB, experimental data with results 
based on the modlfeid V n Driest - denser formula; and his results appear to be 
in good agreement with available experimental data. Tne difference in heat 
flux results duo to MABL a KAHL-X5 is too small .for plotting purposes. The 
differences in the heat flux curves in Figure 1 are partially due to the eddy 
viscosity models in MABL, MABL-KB, and BLIMP programs. The formulas involved 

in the different models to compute eddy viscosity are given in Section II and 
in the Appendix. 

Humber: The computer programs MABL, KAHL-KS, and BLIMP impute a 

Laminar Prandtl number by using the method described in Section IX for a gas 
mixture. BLIMP uses a constant value for the turbulent Prandtl 



number, P+ 


- 0.9; whereas, KABL and KAEL-K3 use the formulation of Gebeci 
(Ref. 17), described in Section II to compute a turbulent Prandtl number. 

A comparison of the heat fluon results, usin.?; ?% ~ constant and Pj_ computed 
according to Gebeci, with experimental data is given in reference (ll) 
which was computed by using KABL. Figure 1 shows a comparison of tne heat- 
flux results partially influenced by different turbulent Prandtl number 
computational models. 

SECTION V 
RPGOMjENDATIOMS 

A comparative study of the three rocket nozzle turbulent boundary 
layer computer programs with regard to heat transfer reveals that all of 
them need some improvements and. modifications in their mathematical models 
and. numerical methods before any one program would be acceptable as a 
standard computing procedure for future theoretical investigations. 

The consideration of the numerical solution methods used in these corn- 
outer programs suggests that the BLIMP with its integral matrix solution 
method will require less nodal points and consequently less computer time 
than the Crank -Nicole on implicit finite difference technique used in both 
the KABL & KAHL-KE programs to achieve compatible results and accuracy. 

The BLIKP program considers multi-component diffusion by using bi- 
furcation approximation (Ref. 26); whereas MABL and MABL-KE use a constant 
Lewis number, usually L=Lx=l, which considers the diffusion coefficient to 
neglecting the multi -component diffusion effect. 


be binary. 



The H:\BL~K3 program introduces a turbulent kinetic energy equation to 
calculate eddy viscosity together with a modified Van Dr i e s t - 0 1 au s e r tur- 
bulent model. The BLIMP progreun. lacks in this area end an addition of this 
improved modal to the cods would be appropriate . ELU-'iP should be improved 
with a turbulent number formula given in reference ( 17 ). 

The computer programs need a wall temperature profile as a boundary 
condition and, at present, this prediction is made from scaled test data or 
industry in-house heat transfer calculations which are not accessible. It 
is recommended that a computer program be extended with an option to pre- 
dict a thrust chamber wall temperature profile by coupling the boundary 
layer analysis with the cooling process. 

Considering different aspects of the computer programs available, 

BLIMP appears to be the more suitable one to be improved, as discussed 
above, and used for future heat transfer analysis of rocket engine thrust 
chambers . 
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APPENDIX A 


TURBULENT BOUNDARY LAYER EQUATIONS: including the turbulent kinetic 

energy equation for the Reynolds stress term describing. the flow 
properties of compressible turbulent boundary layers along the wall of 
a rocket thrust chamber are given below. (For detailed derivation see 
Ref. 35). 

When the transverse- curvature effect is neglected, the compressible 
turbulent boundary layer equations in steady state for two-dimensional 
and axisymmetric flows are written in a curvilinear coordinate system as 
follows: 

Continuity.- feCf * E C?* + S'*') ^ J = O C» 


Momentum : 


Energy: 
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dx 
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A2 


V 

where j=o for two-dimensional flow and j=*l for axismmetric flow. 

Symbols are defined in the nomenclature. 

The flow is assumed to be calorically perfect and obeys the 

equation of state, __ __ • . 

p ^ g-RT (p) 

M 

where the time mean of the correlation between fluctuation density and 
temperature, §' t / is neglected because of its small order of 
magnitude. 

TURBULENT TRANSPORT PROPERTIES 

In formulating the turbulent equations of motion (1-4) it was 

assumed that the turbulent flux terms • £ § ^ ue t0 turbulent mass 

transfer, £ <%v ) ^-turbulent momentum transfer, (^ ( >v) / £) Y"* ~~ 

i-i 

turbulent heat transfer could be related to mean , time averaged 
quantities, through the use of the phenomenological mixing length- 
eddy viscosity concept. Accordingly, the eddy viscosity, £ , eddy 
conductivity, A r and eddy deffusivity, D^, are defined as follows: 



~ C = £ C 

£ ©a /®y) 

C€) 


~C§v/£Y.tf. : 
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- A T O , y) C 2 ' V® A) 
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The turbulent Lewis and Prandtl numbers are defined as 


L t = £ C p 3) t Tt 
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TURBULENT KINETIC ENERGY EQUATION 

To compute the eddy viscosity , £ , defined in Bgn (6) the 
following kinetic energy equation is introduced. 

Considering the continuity and momentum equations for compressible 
turbulent boundary layers, the system of equations due to turbulent 
fluctuations are derived and written below in cartesian tensor notation 


form: 
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Equation (21) is simplified as 
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Adding the three equations for i=j=l, 2, and 3, in E<gn (25), and 


l i * / / / 

defining the turbulent kinetic energy, K q = LI LL -+• V V' -4- A/ nj ^ 
the turbulent energy equation is obtained. 
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where K=K q and the tendency- towards isotropy term was assumed to be 


negligible 
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The first, second, and third terms on the right hand side of equation (26) 
are called the turbulent production, diffusion, and dissipation terms, 
respectively. 

According to Prandtl - Wieghart (Ref. .36) the Reynolds stress, — (.*§0) 
the turbulent kinetic energy, K, and the gradient of time mean velocity, 
0lZ/0yare related as 

-<StfV = k§AK" 2 lffl ( 29 > 

where K is a constant and .A is a dissipation length. Two unknown terms 
in equation (26) are replaced by the following two relations due to Rc&fca 
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(Ref. 15, 37), assuming the compressibility effect, of these terms to be 
small; 

-(§ K 0 + 2 P')v' = cc f A K 1/2 
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Substituting equations (29-31) into equation (26) the final form of 
kinet 

•f dev + 


the turbulent kinetic energy equation is derived, 
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The following boundary conditions are used to solve the equations (1) , 
(5), (IS), (19), (20), and (32) simultaneously; 
at the wall, y=o; 
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At the outer edge of boundary layers, 
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When the dissipation length, .A. , is given, the equations may be solved in 
a closed form. Assuming that the compressibility effect on. the 
dissipation length is negligible (Ref. 38, 39) ,/Vis expressed as a 
polynomial of the distance, y, from the wall. 


A = a Cy/S,) fy + b 0*/8 t ) A- o. d (y/Sj) + e (j 


To avoid a negative value of K in the vicinity of the wall the 
boundary layer is separated into the laminar sublayer close to the wall 
and the wake region. Considering that the effect of molecular viscosity 


term is small compared 


-*• O x ; cnu 

terms (j'l ^ are neglected. Then the 

eddy viscosity, <5 , in the wake region is obtained from equations £ , 2 £?_) <J.S 
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In the region very close to the wall the modified Van Driest 


Model (Ref. 18) based upon the Prandtl mixing length theory is applied; 
that is, the inner eddy viscosity is obtained as 


« = *smi 

(37) 

the mixing length,./, is given 
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and A - Van Driest ’s Damping factor includes the effect of suction or 
mass addition and pressure gradient as shwon by Cebece (Ref. 18). 

The inner eddy viscosity, is used from the wall outward until the 
height at which t 0 --"6Vis reached. From that point on to. the boundary 
layer edge the £ c is used. 



